### My script starts here.  This script will track temperature anomalies through time across all cells in an ASCII
### The user defined threshold will be 28 degrees.  This script will calculate for all cells the longest streak above threshold
### And total days above threshold

# First define directories and install any packages that may be necessary

library('SDMTools')

ASCII.dir = '/home1/99/jc152199/Microclimate Statistical Downscale/Spatial TMax/'
out.dir = '/home1/99/jc152199/Microclimate Statistical Downscale/Spatial TMax/daysabove/'
setwd(ASCII.dir)

# Read in a single ASCII file to act as the template

base.asc = read.asc.gz(paste(ASCII.dir,'20070101_Tmax.asc.gz', sep=""))
base.asc = base.asc * 0

# The following command generates a two-column dataframe of unique row/column positions
# Only at those row/column positions where data is finite (no 'NA' or 'NaN' or infinite values

base.pos = as.data.frame(which(is.finite(base.asc), arr.ind = T))

# Add east and north to this data frame of row/column positions

base.pos$north = getXYcoords(base.asc)$y[base.pos$col]
base.pos$east = getXYcoords(base.asc)$x[base.pos$row]


# Generate a list of ASCII files in the ASCII directory

ASCII.list = list.files(ASCII.dir, pattern='.asc.gz')

# Assign the threshold

threshold=280

# Create a blank data frame to bind data to in a loop                            

clust.data = base.pos[1:1049307,] # Sub-divide base.pos

# Create a for loop to read in ASCII layers
# This loop creates a single dataframe 'clust.data'
# This data frame has a single row for each row/column combination and a single column for each DOY
# The first four rows are X,Y,east, and north,

for (layerx in ASCII.list) 

    {
    
                              
    col.name = paste('t',gsub('_Tmax.asc.gz','',layerx),sep='')
    clust.data[[col.name]] = round(read.asc.gz(paste(ASCII.dir,layerx,sep=""))[cbind(clust.data$row,clust.data$col)]*10)
    cat(layerx,'...\n')
    
    }

#Write out cluster data
    
write.csv(x=clust.data,file=paste(out.dir,'clusterdata1.csv',sep=""),row.names=F)

#The following function was written by JJV
#It will identify numbers in a vector list above a user defined threshhold and group them into unique clusters

extreme.clust = function(x,y,thresh.type){
if (thresh.type=='>=') z = which(x>=y)
if (thresh.type=='>') z = which(x>y)
if (thresh.type=='<=') z = which(x<=y)
if (thresh.type=='<') z = which(x<y)
if (thresh.type=='==') z = which(x==y)
if (length(z)>1){
z.diff = diff(z)
z.clust = NULL; tclust=1
for (i in 1:length(z.diff)){
z.clust = c(z.clust,tclust)
if (z.diff[i]>1) tclust = tclust + 1
if (i == length(z.diff)) z.clust = c(z.clust,tclust)
}
out = rep(NA,length(x))
out[z] = z.clust
return(out)
} else if (length(z)==1) {
out = rep(NA,length(x))
out[z] = 1
return(out)                
} else {
return(NULL)
}
}

threshold = 280 #Set threshold (280 not 28 because all values were multiplied by 10 and rounded to become integer)
pos = 1:365  #Number of days to track pattern across i.e. 1:365
t.data = NULL # Blank data frame to row-bind to

for (i in 1:nrow(clust.data))

    {

    cells.above = extreme.clust(clust.data[i,5:369],threshold,'>=')
    
    cat(i,'.../n')
    
    if(is.null(cells.above))
    
      {
      
      streak.above = 0

      total.above = 0
      
      t.data1 = data.frame(row=clust.data[i,1], col=clust.data[i,2], east=clust.data[i,4],north=clust.data[i,3], streak=streak.above, total=total.above)
      
      t.data = rbind(t.data, t.data1)
      
      #cat(t.data1,'.../n')
      
      }
      
        else
        
        {

        out=aggregate(pos,list(cluster=cells.above),length)

        streak.above = max(out[,2])

        total.above = sum(out[,2])

        t.data1 = data.frame(row=clust.data[i,1], col=clust.data[i,2], east=clust.data[i,4],north=clust.data[i,3], streak=streak.above, total=total.above)

        t.data = rbind(t.data, t.data1)
        
        #cat(t.data1,'...\n')

        }
        
  }
  
write.csv(x=t.data, file=paste(out.dir,'summary6.csv',sep=""),row.names=F)

  
for (ii in 1:nrow(out)) lines( c(out$start.day[ii]-0.5,out$start.day[ii]-0.5+out$length[ii]),c(1,1),col='blue')


